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1.  Introduction 

A  truly  automatic  code  for  ordinary  differential  equations  must  not 
only  handle  the  most  general  case  reasonably  efficiently,  but  must  also 
automatically  detect  those  classes  of  problems  which  are  unreasonably 
expensive  by  general  methods  and  switch  to  methods  which  are  more  efficient 
for  those  problems.  This  paper  will  consider  two  classes  of  problems, 
stiff  and  oscillatory,  for  which  special  methods  can  be  far  more  efficient 
than  general  methods.  This  is  not  to  say  that  there  are  not  other  classes 
of  problems  that  are  worthy  of  special  treatment.  For  example,  linear 
equations  probably  can  be  solved  more  efficiently  if  this  fact  is  known 
[1] ,  but  at  this  time  there  are  no  methods  that  are  sufficiently  more 
efficient  for  linear  problems  than  the  general  methods  that  it  seems  worth 
the  effort  to  detect  linear  problems  automatically.  (Furthermore,  most 
users  can  tell  if  a  problem  is  linear,  while  it  may  be  difficult  for  them 
to  tell  when  it  becomes  stiff  or  oscillatory.) 

Although  it  is  common  to  talk  about  "stiff  differential  equations,"  an 
equation  per  se  is  not  stiff,  a  particular  initial  value  problem  for  that 
equation  may  be  stiff,  in  some  regions,  but  the  sizes  of  these  regions 
depend  on  the  initial  values  and  the  error  tolerance.  For  most  problems 
the  solution  is  initially  in  a  transient  and  an  accurate  solution  demands  a 
stepsize  sufficiently  small  that  the  truncation  error  of  that  transient  is 
small.  For  such  stepsizes  stability  is  not  a  problem.  When  the  transient 
has  decayed  below  the  error  tolerance,  the  problem  may  be  stiff.  At  this 
time  a  stiff  method  must  be  used.  Many  techniques  and  programs  are 
available  for  stiff  equations  ([2],  [3],  [4],  [5],  [6])  so  we  will  not 
repeat  that  material. 

Until  recently  stiff  methods  have  also  been  used  in  the  transient 
region,  but  the  fact  that  they  are  generally  less  efficient  than  nonstiff 
methods  (both  because  of  smaller  error  coefficients  and  the  linear  algebra 
involved)  has  encouraged  several  people  to  investigate  automatic  detection 
of  stiffness. 

The  problem  of  highly  oscillatory  solutions  has  some  parallels  to  the 
stiff  problem.  Again,  the  solution  may  not  be  nearly  periodic  initially, 
but  after  a  transient  starting  phase,  may  tend  towards  a  periodic  solution 
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or  have  a  nearly  periodic  behavior.  There  are  some  methods  that  are 
applicable  in  the  latter  phase,  for  example,  Mace  and  Thomas  [9],  Graff 
[6] ,  Graff  and  Bettis  [7]  and  Petzold  [10].  However,  these  methods  cannot 
be  used  in  the  transient  phase  so  it  is  essential  to  detect  the  time  when 
nearly  periodic  behavior  has  begun  and  to  estimate  the  period  reasonably 
accurately.  There  are,  of  course,  problems  for  which  the  user  knows  that 
the  solution  is  nearly  periodic  throughout.  Satellite  orbits  are  a  case  in 
point.  (Most  of  the  early  methods  were  developed  for  these  problems.)  In 
such  cases  the  period  is  known  reasonably  accurately  so  there  is  no 
detection  problem.  This  paper  is  particularly  concerned  with  those 
problems  which  may  become  nearly  periodic  in  later  stages  and  methods  for 
detecting  this  behavior  in  order  to  switch  to  an  appropriate  scheme. 

Methods  for  nearly  periodic  problems  are  generally  known  as 
mult  Revolutionary  from  their  celestial  orbit  background.  The  idea  of  such 
methods  is  to  calculate,  by  some  conventional  integrator,  the  change  in  the 
solution  over  one  orbit.  If  the  period  of  an  orbit  is  T  (for  a  moment 
assumed  fixed),  then  a  conventional  integrator  is  used  to  compute  the  value 
of 

d(t)  =  y(t  +  T)  -  y(t) 
by  integrating  the  initial  value  problem  y'  =  f(y)  over  one  period  T.  If  we 
consider  the  sequence  of  times  t  =  mT,  m  integral,  we  have  a  sequence  of 
values  y(mT)  which  are  slowly  changing  if  y  is  nearly  periodic.  The 
conventional  integrator  allows  us  to  compute  the  first  differences  d(mT)  of 
this  sequence  at  any  time  mT.  Under  appropriate  "smoothness"  conditions 
(whatever  that  means  for  a  sequence)  we  can  interpolate  or  extrapolate  for 
values  of  d(mT)  from  a  subset  of  all  values  of  d,  for  example  from  d(kqT), 
k  =  1,  2,  3,...,  where  q  is  an  integer  >  1,  and  thus  estimate  y(mT)  by 
integrating  only  over  occasional  orbits. 

In  a  satellite  orbit  problem  it  is  fairly  easy  to  define  the  meaning 
of  "one  period."  For  example,  one  could  use  a  zero  crossing  of  a  particular 
coordinate,  or  even  a  fixed  period  based  on  a  first  order  theory.  In  her 
thesis,  Petzold  considered  problems  for  which  it  is  difficult  to  find 
physical  definitions  of  the  period  and  examined  a  method  for  determining 
the  approximate  period  by  minimizing  a  function  of  the  form 


-  3  - 


I(t,  T)  =  |  |y(T  +  T)  -  y(T)  |  |t 

where  the  norm  measures  the  values  of  y(r  +  T)  -  y(t)  approximately  over 
the  range  Tf:(t,  t  +  T).   The  actual  norm  she  used  was 

I(t,  T)  =  /    y(T  +  T)  -  y(T)  L  dx 
t  z 

where  T  was  the  last  estimate  of  the  period.  The  use  of  T  was  for 
pragmatic  reasons.  Ignoring  that  detail,  the  value  of  T  which  minimizes 
I(t,  T)  is  a  function  of  t,  and  T(t)  was  said  to  be  the  period  of  the 
solution.  This  enabled  d(t)  =  y(t  +  T(t))  -  y(t)  to  be  calculated  and 
multirevolutionary  methods  to  be  used.  The  variable  period  was  handled 
easily  by  a  change  of  independent  variables  to  s  in  which  the  period  is 
constant,  say  1.   The  equation 

t(s  +  1)  -  t(s)  =  T(t) 
was  appended  to  the  system 

z(s  +  1)  -  z(s)  =  g(s) 
where  z(s)  =  y(t(s))  and  g(s)  =  d(t/s)  for  integer  values  of  s.   (When  T  is 
constant,   this   is   the  analog  of   the  old  device  for  converting  a  non- 
autonomous  system  to  an  autonomous  system  by  appending   the  differential 
equation  t'  =  1.) 

The  scheme  for  period  calculation  used  by  Petzold  suffers  from  three 
drawbacks.  The  first  drawback  is  that  it  is  fairly  expensive,  involving  a 
numerical  approximation  to  the  first  two  derivatives  of  I(t,  T)  by 
quadrature  which  requires  integration  for  y(x)  over  two  periods.  In  the 
experimental  implementation,  integration  was  repeated  for  every  iteration 
of  a  Newton  method  to  minimize  I(t,  T)  by  solving  8I/8T  =  0.  This  could 
have  been  eliminated  by  saving  all  values  and  interpolating,  but  the 
storage  cost  becomes  high  and  the  quadrature/interpolation  cost  remains 
non-negligible.  The  second  drawback  is  that  a  reasonably  accurate  period 
estimate  is  needed  for  the  Newton  iteration  to  converge.  Outside  the 
region  of  convergence  for  Newton's  method  a  search  scheme  for  a  minimum 
could  be  used  but  this  would  be  very  expensive  because  of  the  computation 
involved  in  each  quadrature  even  if  all  previously  computed  values  could  be 
saved.    This  makes   the  approach  very  unattractive   for  initial  period 
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detection  when  there  is  no  starting  estimate.  The  third  drawback  is  that 
minimizing  a  function  subject  to  several  sources  of  error  (including 
truncation  errors  in  the  integration  and  quadrature,  and  roundoff  errors 
revealed  by  considerable  cancellation  in  ||  |y(T  +  T)  -  y(f)||)  is  likely  to 
yield  a  fairly  inaccurate  answer.  Since  the  value  of  d(t)  =  g(s)  is  quite 
sensitive  to  small  absolute  changes  in  the  period  T  which  may  be  large 
relative  to  the  period,  the  function  g(s)  may  not  appear  to  be  very  smooth. 

This  paper  discusses  an  alternate  approach  to  the  period 
identification  problem.  It  overcomes  the  cost  and  convergence  problems, 
and  also  seems  to  help  with  the  sensitivity  problem.  This  is  discussed  in 
the  next  section  along  with  stiffness  detection.  The  third  section  reviews 
multirevolutionary  multistep  integrators  and  a  technique  for  handling 
variable  periods  based  on  the  period  identification  algorithm.  The  fourth 
section  discusses  a  numerical  example  while  the  final  section  discusses 
unsolved  problems. 
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2.  Periodic  and  Stiffness  Detection 

A  fully  automatic  method  should  be  able  to  detect  problems  with 
special  properties  that  can  be  solved  more  efficiently,  but  the  cost  of 
detection  should  be  low  compared  to  the  integration  cost  so  that  problems 
without  those  properties  do  not  cost  appreciably  more.  Since  stiffness  and 
nearly  periodic  behavior  are  properties  that  may  appear  at  any  point  in  the 
solution,  the  detection  process  must  operate  continuously.  If  it  is  to 
have  a  low  cost,  it  must  not  take  more  than  a  few  operations  on  available 
intermediate  or  final  results  of  a  standard  integrator.  This  section  first 
discusses  periodic  behavior  detection,  then  stiffness  detection. 

We  have  been  deliberately  imprecise  about  the  meaning  of  "nearly 
periodic,"  and  will  continue  that  way  with  the  working  definition  In  our 
minds  of  "the  type  of  problem  that  can  be  handled  efficiently  by 
multirevolutionary  methods."  The  types  of  problems  that  have  the  required 
properties  are  differential  equations  for  which  there  exist  functions 
F(t,  t)  and  T(t)  such  that 

y(t)  =  F(t,  t) 

is  a  solution  of  the  differential  equation  dy/dt  =  f(y,  t)  ,  F  is  periodic 
in  t  with  period  T(t),  that  is, 

F(t  +  T(t),  t)  =  F(t,  t), 

for  all  t  and  t,  9F/9t  is  very  small  compared  to  3F/8t,  and  T(t)  is  slowly 
varying.  Here,  t  and  t  are  the  "fast"  and  "slow"  times  as  shown  in  Figure 
1.  The  "meaning"  of  this  representation  is  that 
P(t)  =  {F(t,  t),  te[0,  T(t)]}  is  the  local  periodic  behavior  of  the 
solution,  and  the  change  of  P(t)  with  respect  to  t  represents  the  way  this 
behavior  slowly  changes.  (This  representation  is  only  valid  for  problems 
which  are  not  phase  locked  to  a  periodic  driving  function.)  F(0,  t)  was 
called  a  quasi-envelope  by  Petzold.    It   is   the  function  z(t)  defined 
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earlier  for  a  discrete  set  of  points  only 

t 


olution 


T(t) 


periodic 
direction 


Figure  1.   Nearly  Periodic  Solution  Family 

This  representation  is  not  unique,  but  depends  on  the  choice  of  the 
period  T(t)  and  the  values  of  F(0,  t)  over  the  initial  period.  It  is 
convenient  to  consider  a  change  of  variables  to  (s,  t)  with  s  =  x/T(t).  In 
the  new  coordinate  system,  F(s,  t)  has  period  1  in  s  and  a  unique  quasi- 
envelope  is  defined  for  any  fixed  s  in  terms  of  F(0,  t). 

The  "period"  of  a  nearly  periodic  function  has  not  yet  been  defined. 
We  could  use  some  intuitively  reasonable  mathematical  description,  in  which" 
case  we  would  have  to  seek  computational  algorithms  for  its  approximation. 
However,  the  period  is  most  easily  defined  in  terms  of  the  algorithm  used 
to  calculate  it.  It  should,  of  course,  yield  the  exact  period  for  periodic 
functions  and  be  close  for  small  perturbations  of  periodic  functions.  This 
replaces  an  analysis  of  the  accuracy  of  period  calculation  with  an  analysis 
of  the  efficiency  of  the  multirevolutionary  method  with  respect  to 
different  period  definitions.   This  latter  may  be  an  easier  task. 

Petzold's  period  definition,  based  on  minimizing  the  norm  in  eq.  (1), 
is  very  expensive  to  apply  and  cannot  be  considered  as  a  technique  for 
determining  if  an  arbitrary  output  of  an  integrator   is  nearly  periodic. 
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Therefore,  we  look  for  alternate  definitions  of  the  period.  First,  note 
that  if  the  oscillation  is  due  to  a  periodic  driving  function,  we  probably 
know  its  period  or  can  examine  the  system  which  generates  the  driving 
function  directly.  Hence,  we  can  restrict  ourselves  to  autonomous  systems 
or  systems  which  when  made  autonomous  are  nearly  periodic.  (This  means 
that  the  partials  of  the  system  with  respect  to  time  are  small. )  The 
solution  of  an  autonomous  system  is  completely  determined  by  the 
specification  of  the  value  of  the  solution  vector  y_  at  one  time.  That  is 
to  say,  if  we  identify  two  times  on  the  solution  such  that  y_(t^)  =x(t2^' 
we  know  that  the  solution  is  periodic  with  period  to  -  ti.  This  suggests 
determining  the  period  by  looking  for  minimum  of  | y_(t ^ )  ~  7.^2^11*  ^he 
cost  of  this  is  not  particularly  low  and  it  requires  a  clever  adaptive 
program  with  a  lot  of  heuristics.  The  form  of  the  program  is  to  choose  an 
initial  point  t,  and,  as  each  new  value  y^(t)  is  calculated,  to  see  if 
|y_(ti)  -  y_(t)|  has  passed  a  local  minimum.  If  so,  interpolation  is  used 
to  locate  the  minimum.  If  the  value  of  the  norm  at  the  minimum  is  small 
(first  heuristic  parameter),  we  have  possible  periodic  behavior.  If  not, 
we  must  continue  to  examine  more  y(t).  However,  since  the  periodic 
behavior  of  y  may  not  have  started  by  t-.,  we  must  also  advance  t-, 
occasionally  (next  heuristic  parameter).  After  some  experiments,  we 
abandoned  this  approach. 

Another  way  of  defining  the  period  is  to  identify  certain  points  on 
the  solution  at  which  a  simple  characterization  is  repeated,  such  as  zero 
crossing.  The  solution  itself  may  not  have  zero  crossings  and  if  it 
consists  of  a  periodic  function  superimposed  on  a  slowly  growing  function, 
it  may  be  difficult  to  choose  any  value  which  is  crossed  periodically. 
However,   its   derivative  will  have   periodic   sign  changes,   so  we  have 

experimented  with  a  definition  of  period  based  on  the  zero  crossings   of 

T  *■  T 

c_  y_   where  c^   is   the   transpose   of  a  vector  of  constants.   The  program 

T 
examines  the  integrator  output  for  positive-going  zero  crossings   of  _c  y_'. 

(Currently,   c_  is  a  vector  of  the  weights  provided  by  the  user  for  error 

norm  calculations.)  Anything  but  a  simple  periodic   solution  may   lead   to 

more   than   one   zero   crossing   in  a  single  period,   so   the  norm 

|y'(tj_)  -  y'  (t2 )  I   is  also  examined,  where  ti  and  t^   are  a  pair  of   zero 

crossings.    If  the  norm  is   small,   the   possibility  of  a   period   is 
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considered.   The  procedure  used  is  as  follows: 

T 

1.  Identify  a  positive  going  sign  change  in  c y  . 

2.  Interpolate  to  find  the  t  value  t  t  of  the  zero  crossing.  Also 
compute  interpolated  values  of  y_,  y_'. 

3.  Save  these  values.   (Up  to  ten  prior  values  are  saved.) 

4.  Compare  current  value  with  each  prior  value  in  turn  until  a  small 

y'^i  j  -  y'  .„...,„..   is  found. 
i  ii  old  J-   current1  ' 

5.  Save  period,,^,  =  t  _„____«.  -  t^1  A. 

r  new    current    old 

6.  Compare  period  with  period  ■,  .  if  one  has  previously  been 
calculated.  If  they  are  relatively  close,  accept  the  new  period  and 
switch  to  multirevolutionary  methods. 

7.  Examine  several  backward  differences  of  recent  periods.  If  they  seem 
to  be  smoothly  varying,  accept  new  period. 

The  last  test  was  found  to  be  necessary  for  some  problems  with  variable 
periods. 

As  can  be  seen,  there  are  numerous  heuristics,  which  implies  that  much 
tuning  is  possible.  However,  it  is  important  to  note  that  tuning  effects 
efficiency  only.  If  the  tests  for  periodicity  are  too  stringent,  the 
standard  integrator  will  continue  too  long;  if  they  are  too  lenient, 
multirevolutionary  methods  will  be  invoked  before  they  are  efficient. 
However,  since  they  will  then  run  using  a  stepsize  of  one  period,  they  will 
perform  very  little  worse  then  the  conventional  integrator.  (The  only 
losses  are  additional  overhead. ) 

The  multirevolutionary  method,  described  in  the  next  section,  is  a 
modification  of  a  standard  integrator.  It  calls  on  a  subroutine  to 
evaluate  g(z)  =  z(s  +  1)  -  z(s)  given  z(s).  It  can  suffer  from  stiffness 
in  exactly  the  same  way  that  an  ordinary  integrator  can  suffer  from 
stiffness:  if  H9g/9z  is  large  the  method  may  be  unstable  and  the  corrector 
iteration  will  not  converge  unless  a  Jacobian  J  =  3g/8z  is  used  in  a  Newton 
iteration.  Shampine  [11]  has  suggested  monitoring  the  size  of  the  Jacobian 
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by  estimating  its  norm  when  two  or  more  function  evaluations  are  used  in  a 
single  step.   Essentially, 

|  |  g(z-,  )  -  g(z2)| 

L  =  max  — 


lfZl  '  Z2|| 

is  calculated,  where  the  max  is  taken  over  all  steps  and  z^  and  z2  are  two 
different  values  of  z  for  which  g  is  evaluated  in  one  step.  (In  practice, 
it  seems  preferable  to  take  an  exponentially  weighted  max  such  as 

,r,   n  T     I  I  sCz,  )  -  g(z2)|  I 

L    =  max  (0.9  L  -,., * 4 ) 

new       v     old'     z^  -  z2 

but  this  is  yet  another  tuning  heuristic.) 

This  technique  can  be  used  in  both  the  regular  integrator  used  to 
calculate  g  and  the  multirevolutionary  integrator.  However,  there  are  a 
number  of  tuning  questions  that  pose  some  difficulties.  One  could  decide 
to  switch  to  a  stiff  method  the  moment  that  J  becomes  large  enough  to 
restrict  the  stepsize  in  which  case  there  are  no  problems.  However,  that 
is  not  efficient  because  nonstiff  methods  are  both  considerably  lower  in 
cost  per  step  and  have  a  smaller  error  tolerance.  The  natural  solution  to 
consider  is  to  continue  with  the  nonstiff  method  until  the  estimated 
stepsize  that  could  be  used  in  a  stiff  method  is  sufficiently  larger  to 
offset  the  increased  cost  per  step.  This  requires  estimating  the  error  in 
stiff  methods  by  estimating  various  derivatives.  This  is  done  with 
suitable  difference  formulas  but  if  care  is  not  used,  the  derivatives 
estimated  may  appear  to  be  such  that  little  step  increase  is  possible  with 
a  stiff  method.  The  reason  for  this  difficulty  is  that  most  codes  do  not 
directly  restrict  their  stepsize  on  the  basis  of  a  Jacobian  estimate  and 
stability  needs.  At  the  most  they  restrict  h||j|  so  that  the  corrector 
would  converge  if  iterated.  However,  higher  order  methods  may  be  on  the 
verge  of  instability.  Once  h||j|  is  too  large  for  stability,  small  errors 
grow  rapidly.  The  automatic  error  control  quickly  reduces  the  step  to  keep 
the  error  near  the  tolerance  level  and  a  well-engineered  nonstiff  solver 
will  produce  a  perfectly  good  solution,  albeit  slowly.  However,  the 
solution  contains  errors  of  the  order  of  the  error  tolerance  due  to  these 
marginally  stable  components,  and  these  errors  usually  oscillate.    When  a 
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difference  formula  is  applied  to  them,  large  values  result,  and  completely 
obscure  the  derivatives  we  want  to  estimate.  For  example,  the  marginal 
stability  could  introduce  an  error  of  (-l)ne  into  the  numerical  solution  y 

at  the  n-th  step.    If  we  form  the  k-th  backward  difference  to  estimate 

k  (k^  n  k 

h  yv  '   we  have  a  component  due  to  this  error  of  (-1)  2  e.   If  we  now  ask 

k  fk") 
what   stepsize  ah  can  be  used  in  a  stiff  method  whose  error  is  C,  h  yv  '    to 

achieve  an  error  of  t,  we  find  that  a  =  (2  Ci  )   '   independent  of  current 

h.   For  BDF,  C^  =  1/k  so  a  is  always  less  than  one,  falsely  indicating  that 

the  stepsize  cannot  be  increased.   To  avoid  this  difficulty  it  is  necessary 

to  keep  h||j|    small  enough   that   components  with   the  most  negative 

eigenvalues  are  at  least  moderately  damped. 

In  addition  to  estimating  |j||,  it  is  possible  to  estimate  the 
largest  eigenvalue  (or  eigenvalue  pair  if  complex  conjugate)  using 
evaluations  of  g(z)  from  more  than  one  step,  as  suggested  in  Gear  [12]. 
However,  experiments  on  that  technique  indicate  that  the  eigenvalue 
estimates  are  not  too  reliable.  The  reason  for  wanting  to  know  the  largest 
eigenvalue  is  to  know  whether  to  use  a  technique  efficient  for  real 
eigenvalues,  such  as  BDF,  or  a  technique  better  suited  to  eigenvalues  close 
to  the  imaginary  axis,  such  as  BLEND  [13]  or  implicit  Runge-Kutta  [3].  K. 
Stewart*  pointed  out  that  it  is  sufficient  to  wait  until  a  decision  to  use 
stiff  methods  has  been  made. 

At  that  time  a  Jacobian  must  be  calculated,  and  power  methods  can  be  used 
to  determine  the  arguments  of  the  large  eigenvalues.  (This  poses  an 
interesting  question  for  the  numerical  linear  algebraist:  how  to  calculate 
cheaply  the  maximum  argument  of  all  eigenvalues  exceeding  a  certain  size  in 
a  matrix. )  At  the  time  the  Jacobian  is  first  calculated,  it  can  also  be 
checked  for  other  properties  such  as  handedness  and  sparsity  so  that  a 
decision  on  which  linear  equation  scheme  to  use  can  be  made. 


*K. Stewart,  Jet  fropulsion  Lab,  Pasadena,  CA.  Private  communication. 
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3.  Variable  Period  Mult  Revolutionary  Methods 

In  the  original  coordinates  we  have  z(t  +  T)  -  z(t)  ■  d(t).  For  small 
T  this  says  z'(t)  =  d(t)/T.  Hence,  it  is  not  surprising  that  the  numerical 
interpolation  for  z(t)  given  a  technique  for  computing  g(t)  ■  d(t)/T  is 
very  similar  to  a  numerical  integration  technique.  In  the  new  coordinate 
system,  the  basic  structure  of  the  program  is  an  outer  integrator  which 
solves  the  equations 

z(s  +  1)  -  z(s)  -  g(z) 
t(s  +  1)  -  t(s)  =  T(t(s)) 
using  an  outer  stepsize  H.  The  method  varies  the  order  and  stepsize  just 
as  an  ordinary  integrator  does.  See  Petzold  [10]  for  details.  It  calls  a 
subroutine  to  evaluate  g  and  T  given  z  and  t.  This  is  done  by  integrating 
the  underlying  ordinary  differential  equation  y'  =  f(y)  starting  from  y(t) 
=  z,  determining  when  a  period  has  elapsed  and  computing  g(z)  -  y(t  +  T(t)) 
y(t).  Both  methods  for  defining  the  period  discussed  in  the  pevious 
section  have  been  tried.  The  first,  minimizing  | |y(t  +  T(t))  -  y(t)|j,  is 
now  easier  to  implement  because  the  left  end,  y(t),  is  fixed.  The  only 
tuning  difficulty  is  to  ignore  intermediate  minima,  and  we  have  done  this 
by  considering  only  values  of  T  starting  from  0.9  of  the  previous  period 
estimate.  (If  T  changes  more  rapidly  over  H  than  this,  either  H  is  too 
large  or  the  nearly  periodic  assumption  is  questionable.)  The  norm  actually 
used  has  the  form 


'  i "  to  I  A-(*)): 


where  vi^  is  the  j-th  derivative  of  the  i-th  component,  and  the  A  .  are 
weights.  It  appears  to  be  best  to  use  A  .  =0,  s  ^  1,  and  A^.,  =  weight  of 
i-th  component  of  error  weight  vector.  This  allows  for  arbitrary 
nonperiodic  linear  functions  to  be  ignored.  Higher  derivatives  can  be 
Included,  but  knowledge  of  them  is  subject  to  larger  errors  due  to  the 
inner  integration. 


The  second  method,  looking  for  a  zero  crossing,  has  a  difficulty:   the 

T  ~ 
function  c   y_'  will  not  necessarily  be  zero  at  y(t)  ■  z.   (It  can  be  shown 

that  it  will  be  zero  except  for  roundoff  error  for  linear  problems.)   This 
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has  been  overcome  by  choosing  a  vector  c^  separately  for  each  period  such 
that  p(t)  =  c  j_'  =  0  at  the  start  of  the  period.  A  future  zero  of  p(t) 
defines  the  length  of  the  period.  A  unique  c_  is  defined  by  choosing  c_  to 
maximize  p'(t)  at  the  start  of  the  period  subject  to  |cj  ■  1.  This  value 
is  chosen  because  it  minimizes  the  roundoff  error  effects  in  determining  a 
zero  of  p(t).   The  value  of  £,  apart  from  scaling,  is 

£  =  z "  -  ft  i,7,2*  1 

This  requires  a  knowledge  of  y"  at  the  initial  point  y(t).  This  is 
available  because  a  Runge-Kutta  starter  which  computes  the  first  four 
derivatives  is  used  for  the  multistep  inner  integrator  [14].  After  c_  has 
been  calculated,  future  positive  going  zero  crossings  of  p  are  examined,  y_ 
and  its  derivatives  are  calculated  by  interpolation  at  the  zero  crossing 
point  t  +  t,  and  |y_'(t  +  t)  -  y_'(t)|  is  evaluated.  If  it  is  small 
enough,  the  period  is  set  to  t  and  g  is  calculated. 

The  variable  period  multirevolutionary  integrator  is  based  on  a 
modified  Nordsieck  scheme.  Each  component  of  z  is  represented  by  the 
history  vector 

a  =  [z,  Hg,  H2g'/2,  H3g"/6,...,  H^"1^!  ]T 

Petzold  has  shown  that  in  this  representation  the  predictor  has  the  form 

ln,(0)  =  ^-1 
where  A  is  the  Pascal  triangle  matrix  except  for  the  first  row  which  is 

[1,  1,  aL(r),  a2(r),...,  ot^Cr)] 

where  r  =  1/H.   She  also  showed  that  the  corrector  takes  the  form 

£n  =^n,(0)  +lu 

where   m   is   chosen   so   that    a„    "satisfies"    the    relation 

— n 

z(sn  +  1)  -  z(sn)  =  g(sn)  and  %_  is  the  conventional  corrector  vector  except 
in  the  first  component  which  is  a  function  of  r  =  1/H.  Petzold  gives  these 
functions  for  generalized  Adams  methods.  (They  are  polynomials  in  r. )  The 
corresponding  functions  for  BDF  methods  are  inverse  polynomials  in  r.   They 
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are 

Order         First  Coefficient  of  l_ 

1  -1 

2  -1/(3/2  +  r) 

3  -1/(11/6  +  2r  +  r2) 

4  -1/(25/12  +  35r/12  +  5r/2  +  r3) 

5  -1/(137/60  +  15r/4  +  17r2/4  +  3r3  +  r4) 

6  -1/(147/60  +  203r/45  +  49r2/8  +  103r3/12  +  7r4/2  +  r5) 

Petzold  suggests  a  linear  combination  of  the  generalized  Adams  and  BDF 
coefficients,  for  example,  r  •  Adams  +  (1  -  r)  •  BDF  so  that  the  method  has 
the  properties  of  BDF  formula  for  large  H  (r  ->  0),  namely  stiff  stability, 
and  the  property  of  Generalized  Adams  for  r  =  1,  namely  the  outer 
integrator  is  exact.  Since  it  is  not  proposed  to  use  BDF  methods  until 
stiffness  has  set  in  (when  H  is  large),  it  does  not  seem  worth  considering 
this  complication. 
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4.  A  Numerical  Test 

Several  example  problems  have  been  constructed  using  the  Van  der  Pol 
oscillator  to  give  a  nonlinear  oscillation.  Typical  of  these  problems  is 
the  following  system  of  four  equations 


ii  "  Q L 


"l   "  U' 


=  -(u,  -  Uo)  +  2(uo  -  (u,  -  Uo)  )u. 


u3 
u4 


-10"3(u3  -  1) 
10"3  sin  10"3t 


All  initial  values  were  zero  so  uo  =  1  -  e~  and  ua  =  1  -  cos  .OOlt. 
u^  and  u£  are  the  solution  and  first  derivative  of  a  Van  der  Pol  oscillator 
oscillation  about  a  level  U3  and  peak  amplitude  about  2uo.  The  period  is 
about  2tt  for  small  Ui  and  steadily  increases  to  about  7.63  for  u<%  =  1. 


The   matrix  Q  used  was 
-1        111 
1-111 
11-11 
111-1 


1/2 


,-1 


(It  is  idempotent  and  Q  =  Q  .)  All  components  of  y  oscillate  after  an 
initial  period.  The  periodic  detector  located  the  oscillation  at  about  t  = 
156.  The  integration  was  continued  to  t  =  10,000  with  an  average  outer 
step  of  about  28  periods.  At  that  stepsize  the  outer  problem  is  quite 
stiff.  The  oscillatory  behavior  is  initially  close  to  sinusoidal  and 
changes  to  the  steep-edged  behavior  typical  of  the  Van  der  Pol  oscillator 
by  t  =  1000.  A  local  error  tolerance  of  10  in  the  inner  integrator 
required  about  400  inner  steps  per  period  at  first,  increasing  to  about 

1200  at  the  end.   The  outer  integrator  took  50  steps  with  local   tolerance 

_3 
of   10    from  t  »  156  to  t  =  10,000,  using  154  inner  integrations  over  one 

period   including   those   for  occasional  Jacobian  evaluations   of   g  by 

numerical  differencing,   an  average  speed  up  of  ninefold  over  the  standard 
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inner  integrator  which  would  have  used  about   10   steps   for   the  whole 
problem. 

Plots  of  the  phases  of  the  solution  are  shown  in  Figures  2  to  5.  In 
all  of  these  figures  the  vertical  scales  for  the  four  components  have  been 
renormalized  to  put  y^  between  2i  -  1.9  and  21  -  0.1  for  I  =  1  to  4. 
Figure  2  shows  the  first  phase  of  the  integration  prior  to  detection  of  the 
oscillation. 


Figure  2.  Initial  Phase  before  Period  Detected 
(For  extraneous  reasons,  only  one  integration  point  in  ten  has  been  plotted 
here,  hence  the  jagged  curves.)  The  amplitude  of  the  oscillation  at  this 
point  is  0.99.  The  shape  of  the  oscillation  at  t  =  156  is  shown  in  Figure 
3.  This  shape  changes  and  grows  in  amplitude  to  3.02  by  t  =  10,000,  as  is 
shown  in  Figure  4.  Figure  5  shows  the  smooth  "quasi-envelope"  z  found  by 
the  multirevolutionary  integrator.  It  was  generated  using  the  50  outer 
integration  steps  so  the  actual  solution  y  is  found  by  superimposing  the 
oscillatory  behavior  of  the  form  shown  in  Figures  3  and  4  and  the 
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appropriate  t  values  in  Figure  5. 


i  .i 


iiHi 


7.1 


e.  t 


6.  I 


5.  1 


Figure  3.   First  Period  in  Multirevolutionary  Integration 


Figure  4.   Fiftieth  (last)  Period  in  Multirevolutionary  Integration 
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a.  2»m.  4«m.  saaa.  saaa.  laaaa. 

isaa.  sa*a.  ca«a.  7eee.  ^aaa.         iiaee. 


Figure   5.      Quasi-envelope   z 


0.  aaaa.  4000.  aaaa.  saaa.         10000. 

iaaa.  sbbb  caae .  ?aaa.  9000.  11000. 


Figure  6.   Quasi-envelope  Corresponding  to  u  (Qz) 
It  should  be  noted  that  these  are  not  small  oscillations.   For  example,  the 
range  of   z,   in  Figure   5  (bottom  line)  is  -1.13  to  -0.6  (approx.  ).   The 
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oscillation  changes  in  amplitude  from  0.99  to  3.02  (peak  to  peak)  over  the 
interval. 

In  order  to  check  the  accuracy,  the  equivalent  quasi-envelope  for  _u 
was  recovered  from  £  by  the  transsformation  _u  ■  Qz_.      Since  uo  and  ut  are 

not   oscillatory,   G3   and  u^    should   also   be  1  -  exp(-10  t)   and 

1  -  cos(10  t),   respectively.   Since  the  cosine  component  is  neutrally 

stable,  any  integration  errors  will  not  be  damped  in  later  steps.   The 

relative  error  in  the  cosine  component  at  t  =  10,000  was  .005  (.008 
absolute) . 

These  results  were  without  tuning.  There  are  a  number  of 
inefficiencies  in  the  software  that  can  be  removed.  For  example,  we  did 
not  give  the  multirevolutionary  integrator  the  information  gathered  during 
the  period  detection  phase  about  g  and  its  differences.  These  can  be  used 
to  allow  a  high  order  start  in  the  multirevolutionary  integrator.  We 
believe  that  additional  improvements  are  possible. 
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5.  Further  Problems 

There  are  a  number  of  additional  problems  of  concern.  Here  we  will 
discuss  three  problems:  non-autonomous  problems,  detecting  the  end  of 
periodicity,  and  the  multiple  oscillator  problem.  Some  problems  require 
only  simple  extensions;  others,  in  particular  the  multiple  oscillator 
problem,  pose  serious  difficulties. 

There  are  two  cases  of  the  non-autonomous  problem  y"  ■  f(t,  y)  to 
consider:  those  in  which  f(t,  y)  is  a  slowly  changing  function  of  t  and 
those  in  which  the  t-dependence  is  responsible  for  the  oscillation — we  say 
that  the  oscillation  is  driven.  In  the  former  case  we  can  conceptually 
convert  to  an  autonomous  system  by  appending  the  usual  t"  =  1.  This  term 
is  slowly  varying  so  the  solution  of  the  enlarged  problem  remains  "nearly 
periodic. "  In  the  latter  case  we  can  determine  the  period  by  examining  the 
driving  term  (that  is,  the  t-dependent  terms  in  f)  and  continue  to  use  the 
same  method. 

Some  nonlinear  oscilltors  are  such  that  a  variable  in  the  system 
increases  to  a  point  that  the  oscillation  is  quenched.  Then  there  is  a 
period  of  relaxation  until  it  starts  again.  For  such  systems  an  automatic 
program  must  detect  the  end  of  the  oscillation  and  revert  to  a  conventional 
method.  This  is  analogous  to  the  problem  of  detecting  a  derivative 
discontinuity  in  the  solution  of  a  conventional  differential  equation  and 
similar  techniques  can  be  used.  When  the  multi-revolutionary  integrator 
calls  for  an  evaluation  of  the  local  period  and  of  g(z),  the  period 
detection  scheme  will  be  unable  to  find  a  period  anywhere  close  to  prior 
value,  or  even  to  find  one  at  all.  After  it  has  looked  a  modest  distance 
beyond  the  expected  value,  it  should  report  failure  so  that  the 
multirevolutionary  integrator  can  reduce  its  stepsize  to  find  the 
"discontinuity"  where  the  oscillator  begins  to  be  quenched.  When  the 
stepsize  has  been  reduced  to  one  of  only  a  few  periods,  the  software  can 
revert  to  a  conventional  integration  method. 

The  multiple  oscillator  problem  poses  difficulties  unless  there  is  a 
large  gap  in  frequencies  between  the  two  highest  frequencies.  In  that 
case,  the  lower  frequencies  can  be  viewed  as  the  slowly  changing  components 
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of  F(t,  t)  and  possibly  the  method  can  be  used  recursively  for  the  second 
highest  frequency,  and  so  on.  If  the  two  highest  frequencies  w,  and  o>2  are 
of  the  same  order  of  magnitude,  the  behavior  will  be  far  from  nearly 
periodic  unless  u^/wo  ■  p/q  for  small  integers  p  and  q.  In  that  case  there 
is  a  subharmonic  w^/p  of  the  two  frequencies  which  can  be  used  as  the 
period.  If  not,  there  does  not  seem  to  be  much  hope  unless  the  oscillators 
can  be  isolated  and  treated  separately  by  the  techniques  discussed  above. 
Suppose  we  can  visualize  the  system  as  consisting  of  two  oscillators  u'  = 
P(u,  y)  and  v"  ■  q(v,  y)  where  y  is  a  slowly  varying  term,  and  a  slow  part 
described  by  y"  =  f(y,  u,  v).  If  f  is  linear  in  u  and  v,  it  is  sufficient 
to  find  the  behavior  of  the  average  of  u  and  v  and  this  can  be  done  for 
each  separately.  However,  if  f  is  nonlinear  in  u  and  v,  we  must  also  keep 
track  of  the  relative  phases  of  the  oscillations  of  u  and  v  so  that  each 
time  f  is  evaluated  on  a  coarse  mesh,  the  correct  relative  phases  can  be 
used. 
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